#include "simulation.h"

static
void init_rho( Simulation *sim )
{
    CHECK_SIM(sim);
    {
        const double Xc  = sim->Lx/2;
        const double Yc  = sim->Ly/2;
        const double R   = (Xc<Yc?Xc:Yc) / 2;
        const double R2  = R * R;
        double **    rho = ARRAY_OF(sim->rho);
        int i,j;
        field_set(sim->rho,0.1);
        
        for(j=0;j<sim->ny;++j)
        {
            const double y       = j * sim->dy;
            const double delta_y = y - Yc;
            const double delta_y2 = delta_y * delta_y;
            
            for(i=0;i<sim->nx;++i)
            {
                const double x        = i * sim->dx;
                const double delta_x  = x - Xc;
                const double delta_x2 = delta_x * delta_x;
                
                if( delta_y2 + delta_x2 <= R2 )
                    rho[j][i] = 0.2;
            }
        }
    }
}

static void save_rho( const Simulation *sim, int cycle )
{
    CHECK_SIM(sim);
    field_save(sim->rho,"temp/rho",cycle,sim->dx,sim->dy);

}

int main( int argc, char *argv[] )
{
    Simulation  *sim   = sim_create(20,25,2.0,3.0);
    const double alpha = 0.1 * min_of(sim->dx2,sim->dy2);
    int          iter;
    
    init_rho(sim);
    save_rho(sim,0);
    
    for(iter=1;iter<200;++iter)
    {
        sim_compute_laplacian_rho(sim);
        field_muladd(sim->rho,alpha,sim->Lrho);
        save_rho(sim,iter);
    }
    
    
    sim_delete(sim);
    return 0;
}

